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Abstract 

We study a method for solving the homogeneous Bethe-Salpeter equa- 
tion. By introducing a 'fictitious' eigenvalue A the homogeneous Bethe- 
Salpeter equation is interpreted as a linear eigenvalue equation, where the 
bound state mass is treated as an input parameter. Using the improved 
ladder approximation with the constant fermion mass, we extensively 
study the spectrum of the fictitious eigenvalue A for the vector bound 
states and find the discrete spectrum for vanishing bound state mass. We 
also evaluate the bound state masses by tuning appropriate eigenvalues 
A to be unity, and find massless vector bound states for specific values of 
the constant fermion masses. 
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1 Introduction 



Bethe-Salpeter (BS) equations [H, |], U for the fermion-antifermion systems in gauge 
theories are important for studying the properties of positronium, mesons and so on. 
The BS equation in the improved ladder approximation was applied to the bound 
state problem in QCD. Especially, the numerical solution for the massless pion was 
given in refs. 0, |J, and they succeeded to reproduce the property of the chiral 
symmetry in the light quark sector. Moreover, the BS equation for the B meson was 
solved in the heavy quark limit. 

Let us explain the solvability of the homogeneous BS equation for the pion case. 
The reason is essentially found in the Nambu-Goldstone's theorem]?], [|. First of 
all, the pion is known to be massless without solving the BS equation. Moreover, 
thanks to the axial Ward-Takahashi identity, one scalar component, S, of the pion 
BS amplitude is found to be identical with the mass function of the quark. Then, 
the other components are obtained by solving the resultant inhomogeneous equation, 
where the component S is treated as an inhomogeneous term and linear algebraic 
techniques are applicable. 

The reason why the homogeneous BS equation for the B meson can be solved is 
that the equation in the heavy quark limit becomes an eigenvalue equation. The BS 
amplitude and the binding energy of the B meson are corresponding to the eigen- 
vector and the eigenvalue, respectively. Thus, the problem is reduced to solving the 
eigenvalue equation, which can be done by linear algebraic techniques. 

How about other cases? Generally, there are two problems for solving the homo- 
geneous BS equation for the massive bound state: 

i) We need the fermion propagator in solving the BS equation. It is natural to 
take account of the quantum corrections to the fermion propagator. Actually 
this is required by the consistency with the chiral Ward-Takahashi identity]^, 
10| , |TTf in QCD. We have to work in the time-like region for the center-of-mass 



momentum g M (on-shell), while we carry out the Wick rotation to the relative 
momentum p M (off-shell). In a certain gauge choice, the fermion propagator 
in the homogeneous BS equation takes the form iSp 1 (p± q/2) = p' ± ^/2 — 
£(— (p± q/2) 2 ), where q 2 is identical to the bound state mass squared and the 
time component p° is pure imaginary. The momenta flowing along the fermion 
and antifermion lines become complex. Thus, we need the mass function 
on the complex plane. 

ii) The homogeneous BS equation depends on the bound state mass complicatedly. 
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It is not any linear eigenvalue equation so that it seems impossible to obtain the 
bound state mass and its corresponding BS amplitudes directly. Even when we 
adopt the constant mass approximation for the mass function, the complexity 
is not avoided. The homogeneous BS equation includes both the linear and 
quadratic terms of the bound state mass simultaneously. 

One way to avoid the above two problems is using the inhomogeneous BS equation 



in space-like region. [12], |13[ In general, the solution of the inhomogeneous BS equation 
has the perturbative expansion series which begins with the inhomogeneous part. The 
inhomogeneous BS equation may be solved by the iterations: the inhomogeneous part 
is taken as a first trial solution, and is substituted into the inhomogeneous BS equation 
in order to make more accurate solution. Repeating this iteration many times, we 
obtain the solution with appropriate accuracy. 

Another way is taking further approximation in addition to the ladder approx- 
imation. When we take the BS kernel to be the instantaneous (Coulomb-like) in- 
teraction, the homogeneous BS equation reduces to the Salpeter equation. |L4| On 
the other hand, in ref. [|15|] the authors work for the Euclidean total momentum and 
solve the BS equation by expanding the BS amplitude in terms of SO (4) orthogonal 
Tschebyshev polynomials and dropping higher order polynomials. They investigate 
various mesons with their approximation. [|i~6|j 

Now, let us ask the following question: "How do we solve the above two problems 
i) and ii) to obtain solutions of the homogeneous BS equation?" As for the problem 
i), the Schwinger- Dyson (SD) equation on the complex plane was solved to obtain the 
mass function. WJ^ |18f This analytic continuation is also done using the SD equation 
in the integral form. [In appendix [A], following ref. [19|], we briefly review how to 



derive the fermion mass function on the complex plane.] Then we concentrate our 
attention on the problem ii). 

Traditionally, the BS equation is solved by regarding a coupling constant as an 
eigenvalue and the bound state mass as an input. Q In this paper, to solve the ho- 
mogeneous BS equation for a given coupling constant numerically, we introduce a 
'fictitious' eigenvalue A by replacing the BS kernel K with XK. We notice that it 
is corresponding to replacing the coupling a with Xa in the (improved) ladder and 
constant fermion mass approximations. By adjusting the input bound state mass 
so as to set the fictitious eigenvalue A equal to unity, we get masses of the bound 
states. As stated above, it is natural to use the full propagator which is determined 
by the SD equation in the same approximation as that for the BS equation. But here 
we concentrate our attention on the numerical method to solve the homogeneous BS 
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equation, and we adopt the constant fermion mass approximation. 



In many literatures |2(J, by regarding the coupling constant as an eigenvalue, the 
BS equation is solved for the fermion- ant if ermion bound state in the fixed coupling 
case. The behavior of the 'eigenvalue' a, which corresponds to our fictitious eigenvalue 



A, is studied in detail. For applying to QCD we should use the running coupling pi 



22|, but we do not know the behavior of the fictitious eigenvalue A. Then, before 
calculating the bound state mass, we study the behavior of A relating it with the 
norms of BS amplitudes. Moreover, we check the existence of the discrete spectrum: 
our choice of the running coupling[[L3| becomes a constant in the low energy region and 
the low energy behavior of the theory would be similar to that of the strong coupling 
QED. For a type of BS equation in the strong coupling QED, it is shown pp|, |23"[ that 
there is no discrete spectrum. 

After studying the behavior of A, we solve the original BS equation for the vec- 
tor bound state numerically, and obtain the bound state mass. As an instructive 
example to check the validity of our method, we also calculate the spectrum of the 
orthopositronium, which is a vector bound state of electron and positron. 

This paper is organized as follows. In section ^| we review the homogeneous BS 
equation. The BS amplitude is defined and is expanded in terms of the invariant 
amplitudes. The normalization condition for the BS amplitude is given. Section |3| is 
devoted to explanation of the method for solving the homogeneous BS equation. We 
extensively study the behavior of the fictitious eigenvalue A in the running coupling 
case in section In section |47|, as a consistency check we study the positronium 
in the weak coupling QED, and we calculate the bound state masses in the running 
coupling case. The discussions are found in section [|. 



2 Homogeneous Bethe-Salpeter Equation 

In this section we show the basic formulations in solving the homogeneous Bethe- 
Salpeter (HBS) equation for vector bound state. Our formulations are rather trans- 
parent especially for the vector case. 



2.1 BS Amplitude 

Let us express the fermion-antifermion bound state of the vector type (J PC = 1 ) 
as \V(q, e)) where is the momentum of the bound state and e M is the polarization 
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vector, [e-q = 0, e = — 1.] The BS amplitude x(p; q, e) is defined by 

(0|T^(xh%)|V%e)> = e ^ x J—^ e^ r x (p; g, e) , (2.1) 

where and r M are the center-of-mass coordinate and the relative coordinate, re- 
spectively: 

= — , = x» - . (2.2) 

The BS amplitude is bispinor and is defined only on the mass-shell of the bound state 
q 2 = My, where My is the mass of the bound state. We usually perform the Wick 
rotation to the BS amplitude, and the time component of the relative momentum p^ 
becomes pure imaginary, while the quantities q^ and e M are Minkowski vectors and 
real. 

The BS amplitude of the vector bound state is expanded by eight invariant am- 
plitudes x 1 , • • •, x s - 

8 

x(p; q, e) = T i(P'' 9, e ) x l (p; q) ■ (2.3) 
i=i 

We note that the dependence of the polarization vector e M is isolated in the bispinor 
base Tj(p; q, e). The following choice of the bispinor bases T^p; q, e) (i = 1, ■ • • , 8) are 
convenient: 

rw, r 2 = ^,p](p.q), r s = il/,2], r 4 = i[^,#), 

i ( 2 - 4 ) 
r 5 = (e-p), r 6 =/(e-p), r 7 = ^(p-g)(e-p), r 8 = 2 |X, #](e • p), 



where g M = q^/My and [a, 6, c] = a[b,c] + b[c,a] + c[a,b]. The invariant amplitude 
X l {P'i q) is a scalar function inp 2 andp-g.Q So, after the Wick rotation, it is convenient 
to introduce the real variables u and x by 

p ■ q = iM v u, p 2 = —u 2 — x 2 . (2.5) 

The time component of p M is iu and the magnitude of the spatial components of p^ 
is x in the rest frame g M = (My,0). The invariant amplitude is a function in u and 
x \ X % — X l ( u i x )- By the charge conjugation properties 

Cx(-p; i, e) T C _1 = -x(p;g,e) , 

CTiC-pj^ e)^- 1 = -^(pj^e), (2.6) 

We are considering the bound states which have discrete spectrum. The quantity q 2 — My 
labels the degree of the radial excitation if we introduce it as x — x(p 2 iP ' 9) Q 2 )- 
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where C = ry 2 7°, all the invariant amplitudes x ) are found to be even functions 
in u:[I3| 

X \u,x) = X \-u,x) . (2.7) 



2.2 Homogeneous BS Equation 

The HBS equation reads 

T X = K X . (2.8) 
The diagram for the HBS equation in the (improved) ladder approximation is shown 



P+q/2 




k+q/2 




p-v/2 k-m 

Figure 1: The Feynman diagram of the homogeneous BS equation ( |2.8| ) 
in the (improved) ladder approximation. The quark lines in the RHS of 
the equation are truncated. 



in fig. III. The kinetic part T is given by 

T(p; q) = S F \p + q/2) ® S F \p - q/2) 
with the tensor product notation 

(A ® B) x = A X B . 
In this approximation, the BS kernel K is given by 

K (P ' k) ~ ° 2 -(p - kf {p-kf ) 7 ® 7 

In eq. ( |2.8|) we use the inner product rule 



i^x(p; 9) 



(2tt) 4 ? 



K(p,k)x{k;q) . 



(2.9) 



(2.10) 



(2.11) 



(2.12) 



In the improved ladder approximation we use the running coupling, the asymptotic 
form of which is given by 



Q'o 



An 



m^oo In fi 



2 ' 



(2-13) 



where ao = 12n/ (HN c — 2Nf) with iV c and Nf being the number of colors and flavors, 
respectively. 



The HBS equation (|2.8|) determines the eight invariant amplitudes x\ ■ • X 8 U P to 
the overall normalization constant as well as the bound state mass My. Multiplying 
eq. ( |2.8|) by the conjugate bispinor base T^p; q, e) = joTi^p*; q, e)^7o, taking the trace 
of spinor indices and summing over the polarizations, we finally obtain 



y 2 dydv 



K i:j (u,x;v,y) x 3 (v,y) , 



T i:j (u,x) x J (u,x) 
where the summation over the repeated indices is promised and 
Tij(u,x) = \ tT [r i (p;g,e)T(p;g)r j (p;g,e 
Kij(u,x;v,y) = 



(2.14) 



i 1 r — 

dcosOY^ ^tr [T i {p;q,e)K(p,k)T j (k;q,e) 



(2.15) 



We omit the explicit forms of Ty(w, x) and Kij(u, x\ v, y) here because they are some- 
what complicated. Due to the charge conjugation property (|2.6|) of the bispinor base, 
the 8x8 matrices and Kij are real and self-conjugate: 



Tij{u,x)^ 



Tij(u,x) T 



Tji(u,x) 



Kij(u,x;v,yy = Kij(u,x;v,y) rj 



Kji(v,y;u,x) . 



(2.16) 



The complex conjugate and the transpose are taken in the 8x8 space. Using the 
property (j2.7|), we are allowed to restrict the integral region of v to positive, v > 0. 
The BS kernel in eq. ( |2.14|) is replaced such that 

dvK ij (u,x;v,y)x J (v,y) = J dv [Kyfa, x;v,y) + K^u, x; -v,y)]x J '(v,y) . (2.17) 

v>0 

Thus, we can treat all the variables u, x, v, y as positive values. 

In order to solve eq. ( 2.14|) numerically we discretize the momentum spaces (u, x) 
and (v,y). The effective parameterizations for the variables u and x are 



u 



x 



x 



(2.18) 
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and similarly for v and y. These new variables are discretized at Nbs points evenly 
spaced in the intervals 

U,Ve [Xu,Au} , 

X,Ye[\ Xl A x ] . (2.19) 
Then, the momentum integration becomes the summation 

J y 2 dydv ■■■ D V D X £ VY 3 ■ ■ ■ , (2.20) 

v> 

where 



»>0 V ' Y 



D v = (A v - Xu)/(N BS - 1) , 

D x = {A x - \x)/{N BS - 1) . (2.21) 

The BS kernel Kij(u, x; v, y) has an integrable singularity at (u, x) = (v, y) which 
comes from the pole in the propagator of the gauge boson (cf. eq. (|2.11|) ). In order 



to avoid this singularity we regularize it by taking the four-point splitting |12 



Kij(u,x;v,y) -> - [k^u, x; v+, y+) + K^u, x; v + , yj) 

+K ij (u,x;v-,y + )+K ij (u,x;v-.,yJ)] , (2.22) 

where v± = exp(V^ ± Djj/A) and y± = exp(F ± D x /^)- 

Now, we are ready to solve eq. (|2.14j) numerically. The HBS equation ( p.!4j ) 



becomes a finite dimensional linear equation in terms of the invariant BS amplitude 
X % - All we have to do is to obtain simultaneously the bound state mass My and the 
corresponding BS amplitude x- 

2.3 Normalization Condition 

First, we introduce the conjugate BS amplitude x(p; q, e) as 

(V(q,eW(ymx)\0) = -e l9X J^- 4 e+ x(p;q,e) . (2.23) 
The relation between the conjugate BS amplitude and the BS amplitude is given by 

X(p;q,e) =7o[x(p*;g,e)] f 7o • (2.24) 

Now, let us fix the normalization of the BS amplitude x- Usually the bound state 
\V(q,e)) is normalized by 

(V(q',e')\V(q,e)) = 2q (2n) 3 5^5 3 (q - q') . (2.25) 
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This condition fixes the normalization of x via its definition (|2.1| ). The convenient 



form for the normalization of \ is given by the Mandelstam formula 

. ,dT, , 



Here we use the bra-ket notation || 



2g M . (2.26) 



ix\Am = Jj^qtr [x( P ; q)(A^){p; q)) (2.27) 

with the trace taken over the spinor indices and all the other inner degrees of freedom. 
Multiplying eq. fl2.26|) by and using q' M (d/dq^) = M v (d/dM v ) we get 

2Mv = {xl M^ lx} ■ (2 ' 28) 

It is easy to find the normalization condition in component form as 

N r BT ■ 

2M V = -| x 2 dxdu ( x *(-u,x))*-J- X >{u,x) . (2.29) 

ZW 1 J OMy 

3 The Method for Solving the Homogeneous 
Bethe-Salpeter Equation 

In this section we study the method for solving the homogeneous Bethe-Salpeter 
(HBS) equation and explain how to apply it to the numerical calculation. 
Let us start with writing down the HBS equation 



T X = K X - (3.1) 

This HBS equation simultaneously determines the mass spectrum My \ 
• • • of the bound states and the BS amplitudes Xiipi <2S e )> X2G ; q, e), Xz(p] q^ e )i • • • up 
to their normalizations. Therefore, to solve the HBS equation is to find a mass My 
of a certain bound state and its corresponding BS amplitude xip'^q^)- However, it 



is difficult to solve HBS equation (3.1) as it is. For definiteness, we choose the suffix 



n as il4 1} < M { y ] < il4 3) < • • 

Usually, the coupling constant a is regarded as an eigenvalue of HBS equation 
( p.l|) for a given bound state mass,[|J where we use the ladder and the constant 
mass approximations simultaneously. The dependence on the coupling constant is 
isolated in the BS kernel and the coupling factors out. As far as we adopt those 
approximations, this is extended to the running coupling case by regarding the overall 
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constant «o given in eq. Q2.13| ) as an eigenvalue. On the contrary, if we consistently 
use the full propagator determined by the SD equation in the same approximation 
as that for the BS equation, the kinetic part T also depends on the coupling a as 
well as the BS kernel K does. Then, the HBS equation ( [3.1|) cannot be regarded as 
a simple linear eigenvalue equation for the coupling. 

The original idea to regard the HBS equation as a linear eigenvalue equation for 
the coupling and the BS amplitude is easily extended to generic cases. By introducing 
a 'fictitious' eigenvalue A we interpret the HBS equation ( |3.1| ) as a linear eigenvalue 
equation for a given bound state mass My: 

Q) T X = K X , (3.2) 

where the kinetic part T is regarded as a weight function. We calculate the bound 
state mass for a given coupling a numerically from this equation. When we consider 
the mass My as an input parameter, the form of eq. ( |3.2|) guarantees the solvability: 
we can solve eq. (|3.2|) by standard linear algebraic techniques for eigenvalue equations, 
and obtain the solutions of the eigenvalue A. The fictitious eigenvalue is essentially 
the same as the coupling constant, A oc a, when we use both the (improved) lad- 
der and constant fermion mass approximations. The resultant HBS equation (|3.2j ) 
simultaneously describes various systems parameterized by the coupling constant a. 

In what follows we first study the basic notions of the eigenvalue equation ( |3.2| ) 
and then we explain the method for solving the original HBS equation ( [3.1|) . We note 
that we have to check the existence of the discrete spectrum of A. As we will show in 
the next section, there exists the discrete spectrum in the ladder and constant mass 
approximations. Here and henceforth we consider such case. In order to distinguish 
the quantum numbers of the eigenvalue equation ( |3.2| ) from that of the original HBS 
equation we attach primes to them like n' or 1', 2', 3', • • • when we solve eq. (|3.2| ) 
to obtain the spectrum of A for a given My. The obtained eigenvalue A is a function 
of the bound state mass parameter My, so we write the dependence explicitly as 
X n i(My). For the convenience of the numerical calculation we define the numbering 
of {A„/(My)} so as to satisfy 

< Xv(My) < A 2 '(My) < Xy{My) < ■ • • , (3.3) 

where the equality holds when the eigenvalues are degenerate. In our numerical 
calculation we encounter the accidental degeneracy of the two eigenvalues at some 
point My. We call this phenomena 'level crossing'. 

There are two kinds of norms of the BS amplitude in eq. (|3.2|) . The quantity 
(Xn'\dT/dMy \xn')i which we call T'-norm, gives the normalization for the BS am- 
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plitude as in eq. (|2.28 ). The non-positive T'-norm states are ghosts and cannot 



be normalized. So, we consider only the states which have positive T'-norm. The 
quantity {Xn'\(—T)\xn'}, which we call T-norm, gives also a natural norm for the 
eigenvalue equation (|3.2|) in the following sense. In the case where the eigenvalue A 
is not always real, it is easy to derive the well-known relation in linear algebra: 

(I/A* - l/\* m ,)( Xm '\(-T)\xn>} = . (3.4) 

This relation says that the eigenvalue A is real otherwise the corresponding state \ 
has zero T-norm, and the state \xn'} can form a complete orthonormal set under the 
inner product (x,V0 = (x\(-T)\i>j ■ 

For the purpose that we find the solutions with A = 1, it is enough to consider 
A is real and positive. [We do not consider the possibility that A becomes negative, 
because it implies the repulsive interaction which cannot form any bound states.] 
Then, we throw away all non-real eigenvalues and its corresponding eigenvectors. 
Further, we do not consider the complex My in solving eq. (|3.2| ). If My is complex, 
the corresponding state is complex ghost and is physically unacceptable. 

Now, let us consider the relation between the linear eigenvalue equation (|3T2|) and 
the original HBS equation ( [3.1|) . What we obtain from ( |3.2|) is the set {\ n ',Xn'} 
for a given My , but what we have to find out is the set { My ,Xn}- If we find an 
eigenvalue A n > and its corresponding eigenvector Xri such that A n / = 1 for a certain 
value of the mass My, this solution also satisfies the original equation fl3.1|) : 



= (T - \ nl K) Xn > = {T- K) Xn > ■ (3.5) 

Namely, the set {My,x n '\^n' — 1} is nothing but a solution of the HBS equation 
( ft.ip . Of course, the other states m' ^ n' which do not satisfy X m i = 1 are out of 
the case. Then, the physical spectrum {My^ \ n G N } of the bound state masses is 
determined by the following set of intrinsic relations 

\y{My) = 1, \ 2 ,{My) = 1, \y{My) = 1, • • • . (3.6) 

Clearly, all the roots of these equations give the solutions of the original HBS equation 

Next, we explain how to input the mass My in eq. (|3.2j ) by a systematic manner 
for tuning the eigenvalue to unity (A„/ = 1) numerically. For this purpose we first 
consider the differentiability of the eigenvalue A(My). In general the function X n /(M V ) 
defined by eq. ( |3.3| ) is not always differentiable at level crossing point. We modify the 
definition of the numbering (A„/ — > A n ) so that the functional forms of the eigenvalues 
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become always differentiable. In order to make differentiable form X(Mv) for all range 
of My we invoke the eigenvalue equation ( |3.2| ) itself. Differentiating eq. ( |3.2| ) with 
respect to My, we have|| 

IUA 9T Li 



x(m v ) ; ; m % ■' ; , (3.7) 



where we use the HBS equation for the conjugate BS amplitude. For normalizable 
BS amplitudes we put {x\9T / dMy\x) = 2My according to eq. ( |2.28| ). The equation 
( p.7| ) implies that if the BS amplitude x is normalizable and have non- vanishing T- 
norm, the derivative of the corresponding eigenvalue A exists. Let X n (Mv) (n 6 N) 
denote one of such eigenvalues and let Xn denote its corresponding eigenvector. Once 
we define the numbering (or we can say tagging) X n (My) at a suitable point My, 
we can uniquely identify the functional form of X n (My) according to the "differential 
equation" ( |3.7| ). In other words, the differential equation ( |3.7| ) with a suitable initial 
condition uniquely determines the form of X n (My). 

Here, we explain the systematic method to tune the eigenvalue. We first consider 
the case that the level crossing does not occur, and we identify X n (M v ) = X n i(M v ) 
for n — 1, 2, • • -. The relation ( |3.7|) tells us the response of the eigenvalue X n '(My) for 
the bound state mass My. If we update the bound state mass to My + 5 My, then 
the eigenvalue X n i changes to X n i + 5My(dX n '/dM v ) according to the relation (|3~7|). 
In order to find the root of the equation X n '(My) = 1 numerically, we are better to 
use the Newton's method. The bound state mass is updated according to 



My -> My + 5My 

•\(-T)\Xn>) L 1 
\Xn'j V A n ' 



I dM v 



5M y = l ; n %T\ U ! ■ (3-8) 



When X n '(My) is not differentiable at a level crossing point M , X n i and Xn> in eq. ( |3.8| ) 
should be replaced with A n and Xn respectively. The derivative dX n i/dMy has two 
different values at M depending on how to take the limits, My \ M or My f M Q . 
If the difference between these two values is not so large, it is not a problem for the 
purpose of updating X n r. we do not need the precise value of My so as to satisfy 
X n '(My) = 1 at one time, whereas we need only more suitable value of My than the 
previous trial value of My. Therefore, we generally use the updating procedure (|3.8|) . 

It is important to study the property of the spectrum of A in order to know which 
of the states the solution {My, Xn'l^n' — 1} is corresponding to, the ground state 
or the the first excited state or the other states. In this section we assume that the 
solutions have positive T-norm below. [Indeed, our numerical calculations given in 
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the next section satisfies this assumption.] Under this assumption, all eigenvalues 
A n (My) (n = 1,2, •■■) in consideration are differentiable and decreasing functions. 
Thus, all A n /(Mv)'s are also decreasing functions, and each equation A ri '(My) = 1 in 
eq. Q3.6|) has only one root at most. Many models suggest that A(My) is a decreasing 
function (see for example, ref. ||). [Smaller value of My means larger binding energy 
which should be realized by larger coupling to form tight bound state. A is roughly 
the same as the coupling.] 

We notice that there is a possibility that some of the equations in ( |3.6| ) have no 
solution, e.g., \ n i{My) ^ 1, for any real input value My. Those A's which are already 
smaller than unity (A < 1) at My = never have the solution of A (My) = 1. In this 
case the eigenvalue and its corresponding state (X n f ,Xn') has no physical correspon- 
dence.! In other word, the state xw is annihilated in the physical states which are 
determined in terms of eq. (|3.6| ). We call this phenomena the "level annihilation". 
We actually encounter the level annihilation in the running coupling case, which is 
studied in the next section. 

Now, we will establish the identifications of the physical states {My \ Xn} with the 
solutions {My, x m /|A m /(My) = 1}. In the numerical calculation these identifications 
are easy tasks. We calculate the set of eigenvalues {A n /} from My = to a certain 
large value. When we use the constant mass and (improved) ladder approximation, 
the point My = 2m is the threshold for the constituent fermions to be liberated from 
the bound state. So, it is enough to calculate A up to My = 2m. Supposing that 
there are / — 1 level annihilations, we obtain Ay < • • • < ^Q-iy < 1 < A;/ < • • • at 
My = 0. Then, the l-th eigenvalue Xy, which becomes unity first as increasing the 
input parameter My and its corresponding BS amplitude xv 9™ e the ground state. 
The (l+l)-th eigenvalue X(i+iy which becomes unity second and its corresponding BS 
amplitude Xm' give the first excited state, and so on. 

Finally, we consider the identification of A n (My) with A m /(My) under the level 
crossing. When there is no degeneracy between any two eigenvalues of X m i 's, we can 
uniquely identify the form of A n (My) according to the continuity: A n (My) = A n /(My) 
for v n G N. [In a model a curve of A (My) terminates at some point My ||24|| , but 
the continuity of the curve does not spoil.] However, in general cases there may 
occur the level crossing. The level crossing is that the accidental degeneracy of states 
for a specific value of the bound state mass My. There exist two possible kinds of 
level crossing: 1) \ n i and A m / are degenerate while both are still differentiable; 2) X n i 



* If we allow My to take complex value, we may always have roots of eq. (3J3). But, in our 
calculation we do not consider complex bound state mass because it is unphysical. 
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and \ m i are not differential) le at level crossing point. Below, we discuss the latter 
kind of level crossing. The important thing here is that X n (My) is defined by the 
differentiability (|3.7|) while X n i(My) is defined by the ordering ( [3.3|) . For simplicity 
here we consider the case when only two states become degenerate at one time with 
different slopes of A (My). If one level crossing between the first and the second 
smallest A's occurs at My = M V i then X n (M v ) = X n /(My) (n = 1,2, 3, •••) for 
M v > M V1 and A 2 (M y ) = X r (M v ), X 1 {M V ) = A 2 >(M y ) and X n (M v ) = A n /(M V ) 
(n — 3, 4, • • •) for My < My\. This identification is easily understood in the general 
cases, where the differentiability of X n (My) uniquely make us identify its form. As a 
result, the phenomena of the level crossing is that two different iable curves X m (My) 
and X n (My) cross at a point My = My cross . 



4 Numerical Calculation 



In this section, we first study the functional form of X(My) in our numerical method. 
We consider the running coupling case (improved ladder approximation) in the con- 
stant fermion mass approximation. Next, we evaluate the bound state mass deter- 
mined by the original HBS equation eq. (p.l|). This study is motivated by the qq 
quarkonia in QCD. Here and henceforth, we rescale all dimensionful quantities by 
Aq CD otherwise stated. A QCD is defined by the blow-up scale of the one-loop running 
coupling in the improved ladder model. 

Let us consider the important things for the numerical calculations. As discussed 
before, in the improved ladder approximation the asymptotic form of the running 
coupling is given by the one-loop renormalization group equation, but we have no 



idea for the functional form in the low energy region. As in ref. [13] we use the form 
of the running coupling: 



a(n 2 ) = 



9V 

47T 



1 

7 
1 

T F 
1 



(t F -t c ) 2 -{t-tcf 



+ 



2t"p(tp 
-t c ) 



2t\ 



if t F < t 
if t c 
if t < t c 



<t <ti 



(4.1) 



where t = ln/i 2 and a>o = 12^/(11^ — 2Nf) with N c and Nf being the number of 
colors and flavors, respectively. In this paper we fix N c = Nf = 3 and t F = 0.5 and 
t c = 2.0. 

In the numerical calculation we have to take care of covering the supports of both 
-T)\x) (T-norm) and (x\dT/dMy\x} (T'-norm). Two supports of the norms 
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should lie within the momentum cutoffs [A[/,A[/] and [Ax, Ax] given in eq. (|2.19| ), 
otherwise we fail to obtain the precise values of the eigenvalue A and its corresponding 
BS amplitudes as well as the bound state masses. Covering the support of T-norm 
is necessary for finding out the correct value of the fictitious eigenvalue A by our 
updating method given in sect. [| The T'-norm plays the role for determining whether 
the state is a physical state (T" > 0) or a ghost state (T" < 0). The ratio of T- and 
T-norms is essential for checking the error 5My defined in eq. (|3.8| ) to converge the 
correct value of the bound state mass. 

For the evaluation of My \ we need to cover the support of (xi'\{~T)\xi'} if there 

(2) 

is no level annihilation. When we evaluate M v we need to cover both the supports 
of (xi'|(— T)\xv} and (x2'|( — T)\x2'}, and so on. This is because the excited state Xn' 
should be orthogonal to the lower states xv, X2', • ", X(n-iy in the sense of eq. (|3.4j ).|^ 
In the following calculations we carefully cover the T-norm and T'-norm supports of 
the first two states, i.e., xv an d X2 1 ■ We show the typical example of the supports of 
T-norm and T'-norm for the choice m = 1.0 and My = 1.2 in fig. |2|. [The condition 



10% 




10% 




Figure 2: The supports of a) T-norm and b) T'-norm for the choice 
m = 1.0 and My = 1.2 with Nbs = 17 for the state xv ■ The upper 9/10 
of each figures is clipped. Horizontal axes of U and X are parameterized 



by eq. (|2.18| ). Summing up all values on the lattice points give us the 
norm for a) and b). 



Ai' (My) = 1 is satisfied by this parameter choice.] 



' If we calculate the decay constant, we also need to cover the support of T'-norm in a similar 
reason. 
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4.1 The spectrum of A 

First, we show the functional form of the fictitious eigenvalue X n i(My) with a choice 
m = 1.0 in fig. |3|. This figure obviously shows that the smallest three A n /(My)'s 

8[' • . ■ 

Xn(M v ) 

4 ' • : . 

2 . ' : . 



o i • • d 

0.5 1 1.5 2 

Figure 3: The functional forms of the smallest three A„/(My). The 
horizontal line indicates \ n >(My) = 1. We use A^s = 17. 



are monotonically decreasing functions and there are no level crossing and level an- 
nihilation. The eigenvalues \y(My) and A2' (My) are the first and the second to 
become unity respectively when the bound state mass increases from zero. Then, the 
states xv and Xi' are corresponding to the ground state and the first excited state, 
respectively. At the free fermion threshold My = 2m all eigenvalues of A become 
zero. There are free-state solutions which satisfy T\ = (eq. ( |3.2|) with A = 0). At 
My = all three slopes of \ n >(My) become flat smoothly, i.e., d\ n >(My)/dMy — > 
as My — > 0. As in eq. (|3.7|), this means that T'-norm vanishes {x n '\dT / dMy\x n '} = 
consistently with the normalization condition (p. 28 ). Actually, the calculated ratio 
of T'-norm over T-norm vanishes at My = 0. This result implies that the discrete 
spectrum at My = is the set of normal states in the sense of the limit My — > 0. 

Second, we study the cases for smaller fermion masses. We show the spectrum of 
the smallest three X n '(My) with m = 0.5 and 0.2 in fig. |j. When we chose m = 0.5, 
the smallest eigenvalue \y(My) takes the maximum value unity at the vanishing 
bound state mass My = as in fig. ^ (a). This means that m = 0.5 is a critical value 
for the level annihilation. If we chose smaller value of the fermion mass (m < 0.5), 
then the level annihilation of the state xv w iU occur. Accordingly, the mass of the 
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ground state changes discretely from zero to My = 0.73 around m ~ 0.5. When 
we chose m = 0.2 as in fig. |] (b), the smallest eigenvalue A]/ (My) is always smaller 
than unity and never reaches to unity over the range < My < 2m. In other words, 
the equation A]/ (My) = 1 has no solution for < My < 2m. This means that the 
state xv has no correspondence with the physical solutions of eq. fl3.1|) , i.e., the level 
annihilation occurs. Then, X2', Xz'i ' ' ' are corresponding to the ground state, first 
excited state and so on, respectively. If the constant fermion mass becomes much 
smaller than 0.2, the next smallest eigenvalue A2'(My) will become always smaller 
than unity as well as Ay (My). When m — 0, maximal number of states will fail to 
have any physical correspondences. 



2.5 

Xn(My) 2 

1.5 

1 

0.5 


1.2 

0.8 
0.6 
0.4 
0.2 




(a) 



0.2 0.4 0.6 0.8 1 



M v 



(b) 



0.1 0.2 0.3 0.4 

Figure 4: The plots of the smallest three A n /(My). We use constant 
fermion masses m = 0.5 in fig. (a) and m = 0.2 in fig. (b). The horizontal 
lines indicate A n '(My) = 1. We use Nbs = 11. 



Here it is important to notice that figs. ||| and |] show that the slopes of the 
smallest three A n /(My) are negative. Further, we actually observe that all states in 
consideration have positive T- and T'-norms, and this leads, by eq. (|3.7|) , to the fact 
that the corresponding A re '(My) (or A n (My)) are decreasing functions. This result 
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supports our assumption that all normal solutions have positive T-norm, then all 
eigenvalues A n /(My) are decreasing functions. 

Let us reconsider the level annihilation from the other side. Since the largest value 
of \ n i(My) is realized at the vanishing bound state mass My = 0, it is enough to 
calculate the eigenvalues Xy, A2', • • • at My = for various constant fermion masses 
in order to see whether the level annihilation occurs. We show the result in fig. [5[ 



Xn' 



2 











0.2 



0.4 



0.6 



0.8 



m 



1 



Figure 5: The smallest three A n ' against m at the vanishing bound 
state mass My = 0. We use Nbs = 14. At each value of the fermion 
mass m we distinguish the three points by the definition of its ordering 
Ai' ^ A 2' ^ A3'. 



This figure shows that the first level xv ls annihilated at m — 0.5, the second level 
Xi' are annihilated at m — 0.16 and the third level X3' at m = 0.13 and so on. In 
other words, the ground state mass is given as a solution of Xy(My) = 1 for m > 0.5, 
while we have to solve \2\My) = 1 for 0.16 < m < 0.5. [The two level crossings 
with respect to the fermion mass seem to occur around m = 0.05 and 0.64.] Recall 
that we are considering only real value of My. If we allowed the bound state mass 
My to take complex value, then we could always find the solution of \ n i(My) = 1, 
and might have no level annihilation. Since we are interested in physical states, this 
solution is beyond the scope of this paper. 

Next, let us consider the existence of the discrete spectrum, which we assume so 
far for simplicity. It is sufficient to study at the vanishing bound state mass My = 0. 
As is seen from eq. (|4.1| ), the running coupling becomes constant in the low energy 
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region, which means the low energy behavior of the model is similar to that of the 
strong coupling QED. For a type of BS equation in the fixed coupling case, there 



is no discrete spectrum of A (or a A) for the tightly bound state, My = 0p0[. So 



it is important to investigate whether the discrete spectrum exists in the running 
coupling case. Figure [| shows that the there exists discrete spectrum for all range 
of the fermion mass m > 0. Clearly the state xv is one °f ^ ne discrete states: if 
we had continuous states, the differences A 2 ' — \i> and A 3 / — X^ would be of order 
0(1/Nbs)- But our result shows the order-one differences. To confirm the existence 
of the discrete spectrum of the fictitious eigenvalue A, we use the three choices of 
lattice sizes N BS = 11, 14 and 17. We show the values of Ai', A 2 ' and A 3 / for the 
fermion mass m = 1.0 in table [I]. Each \ n i changes its value by 0(1/ Nbs) as the 



Nbs 


11 


14 


17 


Ar 


2.24 


2.21 


2.23 


A 2 ' 


6.42 


6.61 


6.66 


A 3 / 


7.85 


8.05 


8.22 



Table 1: The dependences of the fictitious eigenvalues Xy, A 2 / and A3' 
on the lattice size for My = and m = 1.0. The difference A( n+ iy — X n i 
is not order 0{l/N B s) but order 0(1). 



value of N B s is changed, while the differences X n > — X m i are never of such order. This 
clearly shows the existence of the discrete spectrum {A„/}. 

Finally, we note other choices for the running coupling. If we tune A to an appro- 
priate value, say Ao, we obtain the solution of the HBS equation which describes the 
bound states in any system with different gauge group and matter contents. This is 
realized when we use the improved ladder and constant fermion mass approximations 
simultaneously. When we chose other gauge groups and matter content than SU (3) 
and Nf = 3, the parameter a (N c = Nf = 3) = 4n/9 defined in eq. (pj~l"l) is modified. 
Since ao is just an over overall constant in the right-hand-side of eq. ( |3.2|) , we can 
identify 

X = ^a , (4.2) 
Arc 

and we regard that all figures [|, £| and || shows the spectrum of ao(My) scaled by 
47r/9. Figure | shows that for the case 2.23 < Ao < 6.66 the equation Xy(My) = Ao 
has no solution, i.e., the level annihilation occurs and the ground state mass My is 
given by A 2 /(M^ 1) ) = A . For the case 6.66 < A < 8.22 M ( y ] is given by A 3 /(M y ) = A . 
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[The condition Aq > 2.23 corresponds to Nf > 10 for N c = 3.] 



4.2 Bound state mass 

Using the method explained in section ^, we evaluate the mass of the vector bound 
state. We consider both the fixed and the running coupling cases with constant 
fermion mass. The vector bound state in weak fixed coupling case corresponds 
to the orthopositronium, the system of which is well-known both theoretically and 
experimentally |25[1 . We can explicitly check the validity of our method in this system. 



Next, we proceed to the running coupling case. 

In the weak coupling limit the HBS equation for the orthopositronium is solved, 
and we obtain the famous non-relativistic result :p3|, |26| 



2 



where m is the pole mass of the fermion and antifermion. On the other hand, the 
spectrum of A (or Xa) for a small binding energy in eq. ( p.2|) is given by [pO| 




\ n ,(M v ) = —J v - . (4.4) 



Obviously, imposing A n /(My) = 1 in eq. ( |4.4j ) we obtain the non-relativistic result 
( [Op . The smallest eigenvalue Xy corresponds to the ground state (n — 1), and the 
second smallest eigenvalue A2' corresponds to the first excited state (n = 2), and so 
on. Furthermore, the binding energy and BS amplitude have no dependence on the 
gauge parameter, [p^jj f 



Here, we rescale all the dimensionful quantities so as to satisfy the relation ma = 1 
fixing the coupling constant as a = 1/137 when we solve the HBS equation for the 
positronium. When the fermion mass becomes small, the support of the T'-norm 
begins to shift to the infrared region. This is the reason why we use large value of 
the fermion mass compared with the unit scale. 

We check whether the numerical calculation reproduce the relation (|4.3|) in the 

weak coupling limit. The spectrum ( f4.4|) tells us that the quantity Xy/X n ' should be 

an integer. We adjust the bound state mass My as Xi'(My) = 1: we update My until 

the difference 5 My given in eq. fl3.8|) satisfies |<JMy| < ma 4 . We show the spectrum 

of the eigenvalues X n i with Nbs = H, 14 and 17 in table §. The results show that 

II In the fixed strong coupling and constant mass the positronium does not exist and there are 
only continuum spectrum and no discrete spectrum. |p(if In the case of the vanishing center-of-mass 
momentum — the HBS equation is solved and a strong dependence on the gauge parameter is 
observed. 1 28 
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Nbs 


11 


14 


17 


Xy/Xy 


2.12 


2.10 


2.09 


A1//A3/ 


3.13 


3.16 


3.17 


\%i / X41 


3.37 


3.32 


3.29 


A1//A5/ 


4.36 


4.39 


4.38 



Table 2: The spectrum of eigenvalues A n / for N BS = 11, 14 and 17 
when the bound state mass is adjusted such that the largest eigenvalue 
A]/ is equal to unity. This should be compared with the result in the 
non-relativistic limit, Xy/X n i — n - 

Xy/X n i appears as an integer for the first several n's within about 10% errors. As 
for the quantum number n' larger than five, Xy/X n > suffers large numerical errors. 
This simply comes from the fact that we do not cover the supports of such higher 
radial excited states. We calculate the binding energies of the first three states. The 
binding energy is defined by = 2m- and is given in table [3] with Nbs — 11, 
14 and 17. This result should be compared with that in the non-relativistic limit: 



N BS 


11 


14 


17 


B^ [Ryd] 


0.513 


0.489 


0.479 


B& [Ryd] 


0.121 


0.116 


0.115 


B^ [Ryd] 


0.0544 


0.0515 


0.0505 



Table 3: The biding energies for the ground state (B^>), the first excited 
state (B^) and the second excited state (B^) with N B s = H, 14 andl7. 
The values are in the units of the Rydberg energy, ma 2 /2. 

B (n) [Ryd] = l/(2n 2 ) (= 0.5, 0.125, 0.0556, ...), and shows that our method works 
well. 

Now that, we proceed to the running coupling case. We update the bound state 
mass My until the difference 8My given in eq. fl3"gl ) satisfies \8M V \ < 10" 2 . As 
studied in section |4.1| , the ground state mass is given by the solution of Xy(My) = 1 
for 0.5 < m and by X2'(My) = 1 for 0.16 < in < 0.5. In the latter case, the 
ground and first excited states are given by xi 1 an d X3'> respectively, due to the 
level annihilation. We show the binding energies (B^ = 2m- il4 n) ) of the ground 
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and first excited states for various fermion masses in fig. This figure shows that 




2.5 3 

m 

Figure 6: The binding energies 

and for various constant fermion 
masses. The upper half data are of B^ and the lower half are of B^ 2 \ 



_£?(!) = 2m at m — 0.5, i.e., \v is massless state. At m = 0.4 the state \v an d its 
eigenvalue Xy have no physical correspondence due to the level annihilation, and the 
binding energy of only the ground state is calculated from X 2 '(My) = 1. As a result, 
the binding energy of the ground state changes discretely from 1.0 (= 2m) to 0.25 
around m ~ 0.5. We make some comments here. The binding energies B' 1 ' and B^ 
approach to the values — > 1 and B^ — > 0.4 asymptotically for large fermion 
mass. It is plausible that the binding energies of a heavy bound states do not depend 
on the mass of their heavy constituents thanks to the decoupling theorem ||29||. We 
may also read off the values of fermion masses with which we obtain massless bound 
states from fig. ||. 



5 Discussions 



As discussed in the Introduction, ideally we are better to solve the SD equation in 
the same approximation as the BS equation in order to obtain the full propagator 
(or mass function). For the purpose that we study the numerical method for solving 
the HBS equation for the massive bound state, we use the constant, i.e., tree level, 
mass S(x) = m in this paper. It is true that this approximation is good in the weak 
coupling QED and retains the non-relativistic result. For applying our method to the 
real QCD case, however, we should use the mass function in the fermion propagator. 
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Even when we consider the heavy quarkonia, the quantum corrections to the mass 
function have the same order as the binding energy. Especially in the case of the light 
quarkonia (p meson) the effect of the mass function becomes more important. As 
a result, the running effect of mass function is always important for any quarkonia, 
and is same order of magnitude as the binding force. 

In the real QCD there exists the states with negative binding energy (for highly 
excited states) because of the confining potential. Our choice of the running coupling 
does not generate the confining force, while a suitable form of the running coupling 
generates it. The non-relativistic quark model with the Richardson potential gives 
a good result for the low-lying spectrum of the heavy quarkonia. |$0| In the constant 
mass approximation the qq free quark threshold opens at q 2 = (2m) 2 .PI Thus the 
states with negative binding energy are unstable. We observe that infinitely many 
number (the order of total lattice points) of states are lying in the region My < 2m. 
We consider only the states which have positive binding energy. Only a first few state 
will be relevant in the constant mass approximation. 

In this paper, we study the the solutions of the homogeneous Bethe-Salpeter 
equation in the constant mass and improved ladder approximations. We leave the 
application of our method to the quarkonia for future work. 



Appendix 

A Schwinger-Dyson Equation 

In this appendix, we briefly review how to derive the fermion mass function in the 



complex plane. Here, following ref. |19fl , we perform the analytic continuation to the 
SD equation in the integral form. 

In order to solve the Bethe-Salpeter (BS) equation we need a fermion propaga- 
tor Sf{p), which is determined by the Schwinger-Dyson (SD) equation. When the 
ladder approximation is adopted for the BS equation, it is natural to use the same 
approximation for the SD equation. In the strong interaction the mass function of 
the quark propagator receives large quantum correction of the same order Aq CD as 
that of binding force. 

** If we take into account of the running effect of the mass function by solving the SD equation, 
the resultant quark propagator has no pole for the real momentum. [TL7L 181 
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The SD equation in the (improved) ladder approximation is given by 



S F \p) =i-m- J^- C 2 g 2 (p, L u - ^ zS F (k) 7 



(A.l) 



(2ir)H 

where = (p — k)^, C% = (N 2 — 1)/ (2N C ) is a second Casimir invariant and m is the 
bare mass of the fermion. The Landau gauge is adopted for the propagator of the 
gauge boson. When we consider QCD, we use the running coupling in the SD equation 
( |A.1|) . Here for defmiteness we adopt the Higashijima-Miransky approximation |2l|, j22| 
to the running coupling, i.e., g 2 (p, k) = g 2 (max(— p 2 , —k 2 )). In considering QED we 
use fixed coupling g 2 (p, k) = e 2 and define Ci = 1. 

In general the fermion propagator is expanded by two scalar functions A and B: 



iS F \p) = A{-p 2 



The mass function E(z) is defined by E(z) = 
the SD equation ( |A.l ) we find A(z) = 1 and 



A{z)/B{z) 



(A.2) 

Substituting eq. ( |A.2|) into 



Vf \ M 3 ° 2 f 

E(z) = m + — - \ 
lbir z Jo 



ydy 



g 2 (max(z,y)) E(j/) 



(A.3) 



max(z,y) y + S 2 (y) ' 

Now, let us discuss the problem i) stated in sec. |1|. When we solve the homogeneous 
BS equation for the massive bound state, we have to use the mass function on the 
complex plane. The mass function £(z) is needed along the set of parabolic curves 
z = (u =F iMy/2) 2 + x 2 , where My is the mass of the bound state. We can restrict 
ourselves to the case Im z > because of H(z)* = £(z*). It is not difficult to 
calculate such a mass function. In refs. flTT| , [TH | the SD equation ( |A.3[ ) is converted 
into the differential equation for carrying out the analytic continuation of the mass 
function. On the other hand, we can easily solve the SD equation (|A.3|) as the integral 
equation itself by the iteration. [B5| After carrying out the analytic continuation, the 
SD equation becomes 

dy- 



Em = m + 



3C, 



167T 2 



+ 



yg\z) 



, yg 2 ( z ) 



dyg 2 (y) 



C{z,oo) 



y + mv) 



(A.4) 



where C(a, b) is the contour from the point a to b and is lying on the parabolic curve 
with Im z > 0. To say more completely, we should take account of the singularity 
of the running coupling. We can solve the SD equation ( |A.4p by the iteration on the 
parabolic curve as usual. Then, we obtain the mass function on the complex plane 
from eq. ( |A.4 ). 
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